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Abstract 

We consider the discrete Couzin-Vicsek algorithm (CVA) [U [191 [36] , which 
describes the interactions of individuals among animal societies such as fish schools. 
In this article, we propose a kinetic (mean-field) version of the CVA model and 
provide its formal macroscopic limit. The final macroscopic model involves a con- 
servation equation for the density of the individuals and a non conservative equation 
for the director of the mean velocity and is proved to be hyperbolic. The derivation 
is based on the introduction of a non-conventional concept of a collisional invariant 
of a collision operator. 
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1 Introduction 



The discrete Couzin-Vicsek algorithm (CVA) [TJ [T9l [36] has been proposed as a model 
for the interactions of individuals among animal societies such as fish schools. The in- 
dividuals move with a velocity of constant magnitude. The CVA model describes in a 
discrete way the time evolution of the positions of the individuals and of their velocity 
angles measured from a reference direction. At each time step, the angle is updated to a 
new value given by the director of the average velocity of the neighbouring particles, with 
addition of noise. The positions are updated by adding the distance travelled during the 
time step by the fish in the direction speficied by its velocity angle. 

For the modeling of large fish schools which can reach up to several milion individuals, 
it may be more efficient to look for continuum like models, which describe the fish society 
by macroscopic variables (e.g. mean density, mean velocity and so on). Several such 
phenomenological models exist (see e.g. [251 EH [35] ) . Several attemps to derive continuum 
models from the CVA model are also reported in the literature [231 [291 EQ], but the 
derivation and the mathematical 'qualities' of the resulting models have not been fully 
analyzed yet. One can also refer to [322 EE] for related models. An alternate model, the 
Persistent Turning Walker model, has been proposed in [18] on the basis of experimental 
measurements. Its large-scale dynamics is studied in [12]. Additional references on swarm 
aggregation and fish schooling can be found in [7] . Among other types of animal societies, 
ants have been the subject of numerous studies and the reader can refer (among other 
references) to [211 S3] , and references therein. 

In this work, we propose a derivation of a continuum model from a kinetic version 
of the CVA algorithm. For that purpose, we first rephrase the CVA model as a time 
continuous dynamical system (see section [2]). Then, we pass to a mean-field version of 
this dynamical system (section [3]). This mean field model consists of a kinetic equation of 
Fokker-Planck type with a force term resulting from the alignement interactions between 
the particles. More precisely, the mean-field model is written: 

e{d t f + u ■ V x f £ ) = -V u ■ (F £ f £ ) + dAuf + 0(e 2 ), (1.1) 
F £ (x, u, t) = v (Id - u ® uj)tt £ (x, t), (1.2) 

Q%x,t)= and f{x,t)= [ vf%x,v,t)dv. (1.3) 

Here f £ (x,u,t) is the particle distribution function depending on the space variable x G 
M 3 , the velocity direction u e § 2 and the time t. d is a scaled diffusion constant and 
Fq{x,u,€) is the mean-field interaction force between the particles which depends on an 
interaction frequency v. This force tends to align the particles to the direction Q £ which 
is the director of the particle flux j £ . the operators V^- and A w are respectively the 
gradient and the Laplace-Beltrami operators on the sphere. The matrix (Id — u <8> w) is 
the projection matrix onto the normal plane to a;, e <C 1 is a small parameter measuring 
the ratio of the microscopic length scale (the distance travelled between two interactions) 
to the size of the observation domain. Here, the relevant scaling is a hydrodynamic scaling, 



2 



which means that e also equals the ratio of the microscopic time scale (the mean time 
between interactions) to the macroscopic observation time. 

The 'hydro dynamic limit' e — > provides the large-scale dynamics of the CVA model 
(in its mean-field version fll.ip - fll.3l) ). The goal of this paper is to (formally) investigate 
this limit. More precisely, the main result of this paper is the following theorem, which 
is proved in section UJ 

Theorem 1.1 (formal) The limit e — > of f £ is given by f° = pMn where p = p(x,t) > 
is the total mass of f° and Q = Q(x, t) G S> 2 is the director of its flux: 



and Mq is a given function of uo ■ Vt only depending on v and d which will be specified 
later on (see l[4.16\) ). Furthermore, p(x,t) and Q(x,t) satisfy the following system of first 
order partial differential equations: 



where the convection speeds c\, C2 and the interaction constant X will be specified in the 



Hydrodynamic limits have first been developed in the framework of the Boltzmann 
theory of rarefied gases. The reader can refer to [HI EED, ED] for recent viewpoints as well 
as to [6l [El [37] for major landmarks in its mathematical theory. Hydrodynamic limits 
have been recently investigated in traffic flow modeling [U [20] as well as in supply chain 
research [31 [H] . 

From the viewpoint of hydrodynamic limits, the originality of theorem 11.11 lies in the 
fact that the collision operator (i.e. the right-hand side of f ll.ll) ) has a three dimensional 
manifold of equilibria (parametrized by the density p and the velocity director Q) but has 
only a one-dimensional set of collisional invariants (corresponding to mass conservation). 
Indeed, the interaction does not conserve momentum and one should not expect any 
collisional invariant related to that conservation. The problem is solved by introducing a 
broader class of collisional invariants, such that their integral (with respect to uj) against 
the collision operator cancels only when the collision operator is applied to a subclass of 
functions. Here, a generalized class of collision invariants is associated with each direction 
Q on the sphere and the corresponding subclass of functions have their flux in the direction 
of Q. We show that such generalized collision invariants exist and that they lead to (\1.7\i . 
In section 14.4} we show that this system is hyperbolic. The detailed qualitative study of 
the system as well as numerical simulations will be the subject of future work. A summary 
of this work can be found in |13j . 




(1.4) 



(1.5) 



d t p + V x ■ (cipfi) = 0. 

p (d t tt + c 2 (tt- V)fi) + A {Id- O <g> tt)V x p = 0, 



(1.6) 
(1.7) 



course of the paper (see \4-41\ ) an d H-63j) )- 
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An important consequence of this result is that the large-scale dynamics of the CVA 
model does not present any phase transition, in contrast with the observations of |36j . 
Indeed, the equilibrium is unique (for given density and velocity director). Therefore, 
the model cannot exhibit any bi-stable behavior where shifts between two competing 
equilibria would trigger abrupt phase transitions, like in rod-like polymers (see e.g. [24] 
and references therein). Instead, the equilibrium gradually shifts from a collective one 
where all particles point in the same direction to an isotropic one as the diffusion constant 
d increases from to infinity. Additionally, the hyperbolicity of the model does not allow 
lines of faults across unstable elliptic regions, like in the case of multi-phase mixtures or 
phase transitions in fluids or solids (see e.g. the review in [22j and references therein). 

With these considerations in mind, a phenomenon qualitatively resembling a phase 
transition could occur if the coefficients ci, C2 and A have sharp variations in some small 
range of values of the diffusion coefficient d. In this case, the model could undergo a rapid 
change of its qualitative features which would be reminiscent of a phase transition. One 
of our future goals is to verify or discard this possibility by numerically computing these 
constants. 

There are many questions which are left open. For instance, one question is about the 
possible influence of a limited range of vision in the backwards direction. In this case, the 
asymetry of the observation will bring more terms in the limit model. Similarly, one could 
argue that the angular diffusion should produce some spatial dissipation. Indeed, such 
dissipation phenomena are likely to occur if we retain the first order correction in the series 
expansion in terms of the small parameter e (the so-called Hilbert or Chapman- Enskog 
expansions, see e.g. [TT]). A deeper analysis is needed in order to determine the precise 
form of these diffusion terms. Another question concerns the possibility of retaining some 
of the non-local effects in the macroscopic model. It is likely that the absence of phase 
transition in the present model is related to the fact that the large-scale limit cancels 
most of the non-local effects (at least at leading order). The question whether retaining 
some non-locality effects in the macroscopic limit would allow the appearence of phase 
transitions at large scales would indeed reconcile the analytical result with the numerical 
observations. A result in this direction obtained with methods from matrix recursions 
can be found in [TDj . Finally, the alignement interaction is only one of the aspects of the 
Couzin model, which also involves repulsion at short scales and attraction at large scales. 
The incorporation of these effects would allow to build a complete continuum model which 
would account for all the important features of this kind of social interaction. 

2 A time continuous version of the discrete Couzin- 
Vicsek algorithm 

The Couzin- Vicsek algorithm considers N point particles in M 3 labeled by k 6 {1, . . . N} 
with positions at the discrete times t n = nAt. The magnitude of the velocity is the 
same for all particles and is constant in time denoted by c > 0. The velocity vector is 
written cwjf where belongs to the unit sphere S 2 = {uj s.t. \uj\ 2 = 1} of M 3 . 
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The Couzin-Vicsek algorithm is a time-discrete algorithm that updates the velocities 
and positions of the particles at every time step At according to the following rules. 

(i) The particle position of the fc-th particle at time n is evolved according to: 

X™ +1 =X% + cu%At. (2.1) 

(ii) The velocity director of the k-th particle, is changed to the director u% of the 
average velocity of the neighboring particles with addition of noise. This algorithm tries 
to mimic the behaviour of some animal species like fish, which tend to align with their 
neighbors. Noise accounts for the inaccuracies of the animal perception and cognitive 
systems. The neighborhood of the k-th particle is the ball centered at X% with radius 
R > and u>j? is given by: 



J k 

I T n \ 
\ J k I 



Jk= £ w i- ( 2 - 2 ) 

j,\X?-X£\<R 



In the Couzin-Vicsek algorithm, the space is 2-dimensional and the orientations are vectors 
belonging to the unit sphere S 1 in R 2 . One can write u% = e * with 9% defined modulo 
2tt, and similarly u)j? = e *=. In the original version of the algorithm, a uniform noise in 
a small interval of angles [—a, a] is added, where a is a measure of the intensity of the 
noise. This leads to the following update for the phases: 

= 01 + & (2-3) 

where 9% are independent identically distributed random variables with uniform distri- 
bution in [—a, a]. Then, uj% + = e lB k . In (36], Vicsek et al analyze the dynamics of 
this algorithm and experimentally demonstrate the existence of a threshold value a*. For 
a < a*, a coherent dynamics appears after some time where all the particles are nearly 
aligned. On the other hand, if a > a*, disorder prevails at all times. 

Here, we consider a three dimensional version of the Couzin-Vicsek algorithm, of which 
the two-dimensional original version is a particular case. Of course, formula (12.21) for the 
average remains the same in any dimension. For simplicity, we also consider a Gaussian 
noise rather than a uniformly distributed noise as in the original version of the algorithm. 
Therefore, our algorithm updates the velocity directors according to: 

< +1 = «Z, (2-4) 

where are random variables on the sphere centered at uj^ with Gaussian distributions 
of variance V2DAt where D is a supposed given coefficient. If the Gaussian noise is 
discarded, the evolution of the orientations is given by 

u n k +1 = Ql (2.5) 

where is the average defined at (12. 2p . 
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Now, we would like to take the limit At — > and find a time-continuous dynamics. To 
do so, we first consider the deterministic algorithm fl2.ll) . (12.51) and following [29], make 
some elementary remarks. First, because \uj^\ = + |, we have (u k +1 — u} k )(u> k +1 + u> k ) = 
0. Therefore, defining u k +1/2 = (u% +1 + u%)/2 and using (12.51) . we have the obvious 
relation: 



, ,n+l n i 

U k ~ U k 1 



At At 



(Id-^ +1/2 ®^ +1/2 )(^-^), (2-6) 



where Id denotes the Identity matrix and the symbol ® denotes the tensor product of 
vectors. The matrix Id — oo k +1 ^ 2 ® oo k +1 ^ 2 is the orthogonal projector onto the plane 
orthogonal to u k +1 ^ 2 . Relation (12.61) simply expresses that the vector uJ k +1 — u k belongs 
to that plane. 

Now, we let At — > 0. Then, the positions X^it) and the orientations uj k {t) become 
continuous functions of time. If we let At — > in (12.61) . the left-hand side obviously tends 
to duJk/dt. The right hand side, however, does not seem to have an obvious limit. This 
is due to an improper choice of time scale. Indeed, if we run the clock twice as fast, the 
particles will interact twice as frequently. In the limit At — > 0, the number of interactions 
per unit of time is infinite and we should not expect to find anything interesting if we do 
not rescale the time. In order to have the proper time-scale for the model, we need to 
replace the tick of the clock At by a typical interaction frequency v of the particles under 
consideration. For instance, in the case of fish, u~ l is the typical time-interval between 
two successive changes in the fish trajectory to accomodate the presence of other fish in 
the neighbourhood. Therefore, we start from a discrete algorithm defined by 



, ,n+l , ,n 

U k ~~ U k HA n+1/2 „ n+l/2w_ n n \ /« v\ 

— =v(Id-w k ®UJ k >){uj k -uj k ), (2.7) 

together with (I2.ip and in the limit At — > 0, we find the following continuous dynamical 
system: 

cuj k , (2.8) 



dt 
duj k 
~dT 



v (Id -uj k ® u k )u k , (2.9) 



where we have used that (Id — u) k ®u k )u k = 0. If the Gaussian noise is retained, then, the 
limit At — * of the discrete algorithm is the following Stochastic Differential Equation: 

dX k 

-^ = cco k , (2.10) 
duj k = (Id — uo k ® uJ k )(uLJ k dt + \/2D dB t ), (2-11) 

where dB t is a Brownian motion with intensity y/2D. Of course, this At — » limit is 
formal but the convergence proof is outside the scope of the present paper. 
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We slightly generalize this model by assuming that v may depend on the angle between 
ojk and u>k, namely v = u(cos9k), with cos^. = uJk ■ UJk- Indeed, it is legitimate to think 
that the ability to turn is dependent on the target direction. If we are considering fish, 
the ability to turn a large angle is likely to be reduced compared to small angles. We will 
assume that v(cos9) is a smooth and bounded function of cos#. 

3 Mean-field model of the discrete Couzin-Vicsek al- 
gorithm 

We now consider the limit of a large number of particles N —>■ oo. We first consider 
the case without Gaussian noise. For this derivation, we proceed e.g. like in [32J. We 
introduce the so-called empirical distribution f N (x,u,t) defined by: 

1 - 

f N (x, ^,t) = -J2 5 ( x - Xk ^ 5 ^ 

k=l 

Here, the distribution w 6 § 2 -> 8{uj,uj') is defined by duality against a smooth function 
(p by the relation: 

(5(u,u'),ip(u)) = ip{u'). 

We note that 8{uj,uj') ^ 8{uj — uj') because the sphere § 2 is not left invariant by the 
subtraction operation. However, we have similar properties of 8 such as 6(u, uj') = 8{uj' , uj) 
where this relation is interpreted as concerning a distribution on the product § 2 x §> 2 . 
Then, it is an easy matter to see that / satisfies the following kinetic equation 

8 t f N + ao ■ V x f N + V. • (F N f N ) = 0, (3.2) 

where F (x, u,t) is an interaction force defined by: 

F N (x, uj, t) = u(cos 9 N ) (Id - uj ® uj)lu n } (3.3) 

with cos^ = uj ■ lj n and uj N (x,uj,t) is the average orientation around x, given by: 

1 1 ' )l j,\X?-x\<R 

If, by any chance, J N is equal to zero, we decide to assign to uj n (x, uj, t) the value uj (which 
is the only way by which uj (x, uj, t) can depend on uj). In the sequel, this convention will 
not be recalled but will be marked by showing the dependence of uj upon uj 

We recall the expressions of the gradient and divergence operator on the sphere. Let 
x = (xi,X2,Xs) be a cartesian coordinate system associated with an orthonormal basis 
(ei,e 2 ,e 3 ) and let (9,<f>) be a spherical coordinate system associated with this basis, i.e. 
X\ = sin 6 cos 0, x 2 = sin 6 sin <fi, x 3 = cos 9. Let also (eg,e^) be the local basis associated 
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with the spherical coordinate system ; the vectors ee and have the following coordinates 
in the cartesian basis: ee = (cos 9 cos 4>, cos 9 sin 0, — sin 9), = (— sin0, cos0, 0). Let 
f(u) be a scalar function and A = A e eg + A^e^ be a tangent vector field. Then: 

VoJ = d e f e e + -A- d+f e$, V w • A = -^—d e (A e sin 9) + --!— <VV 
sin 6* smt/ sm.9 

If the cartesian coordinate system is such that e 3 = uj n , then 

F w = -i/(cos0) sin# e e . (3.5) 
Back to system (13. 2fl - (13.41) . we note that relation (13 .4p can be written 

u N (x,u,t) = r^S^T , J N (x,t) = [ vf N (y,v,t)dydv. (3.6) 



|y-a;|<i?,weS 2 



We will slightly generalize this formula and consider LU N (x,ui,t) defined by the following 
relation: 

.N / ,, ,\ J N (x,t) v / ., f.V, 



u>»(x,u;,t)= ) ' J»(x,t)= K(\x-y\)vf"(y,v,t)dydv , (3.7) 

where i^(|x|) is the 'observation kernel' around each particle. Typically, in formula (13. 6ft . 
i^(|x|) is the indicator function of the ball centered at the origin and of radius R but 
we can imagine more general kernels modeling the fact that the influence of the particles 
fades away with distance. We will assume that this function is smooth, bounded and 
tends to zero at infinity. 

Clearly, the formal mean-field limit of the particle system modeled by the kinetic 
system (13.21) . (13. 3ft . ( 13. 7ft is given by the following system: 

dj + cu ■ VJ + V w ■ (Ff) = 0, (3.8) 
F(x, uj, t) = z/(cos 9) (Id — uj ® u)u)(x, u, t), (3.9) 

J(x,t) f 
v(x,u,t) = - , J(x,t)= K(\x -y\)v f(y,v,t)dydv , (3.10) 

with cos 9 = uj ■ Uj. It is an open problem to rigorously show that this convergence holds. 
For interacting particle system, a typical result is as follows (see e.g. [32]). Suppose that 
the empirical measure at time t = converges in the weak star topology of bounded 
measures towards a smooth function fj(x,u>). Then, f N (x,u>,t) converges to the solution 
/ of (13. 81) - 03. 101) with initial datum fj, in the topology of continous functions of time 
on [0, T] (for arbitrary T > 0) with values in the space of bounded measures endowed 
with the weak star topology. We will admit that such a result is true (may be with some 
modifid functinal setting) and leave a rigorous convergence proof to future work. 

We will also admit that the mean-field limit of the stochastic particle system ( 12 . 1 Of) . 
( 12.111) consists of the following Kolmogorov-Fokker-Planck equation 

dj + cu- VJ + V w ■ (Ff) = DAJ\ (3.11) 
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again coupled with ( 13. 9ft . (13. lQf) for the definition of F and u;, and where A w denotes the 
Laplace-Belltrami operator on the sphere: 



/ 



1 



d e (sm6dgf) + 



1 



sin 9 



sin 2 9 



4 Hydrodynamic limit of the Mean-field Couzin-Vicsek 
model 

4.1 Scaling 

We are interested in the large time and space dynamics of the mean-field Fokker-Planck 
equation (131TD . coupled with (1531) . (13~T0|) . 

So far, the various quantities appearing in the system have physical dimensions. We 
first introduce the characteristic physical units associated with the problem and scale 
the system to dimensionless variables. Let vq the typical interaction frequency scale. 
This means that v{cos9) = vqv'{cos9) with z/(cos#) = 0(1) in most of the range of 
cos 9. We now introduce related time and space scales t and x as follows: t = u^ 1 
and Xq = ct = c/u . This choice means that the time unit is the mean time between 
interactions and the space unit is the mean distance traveled by the particles between 
interactions. We introduce the dimensionless diffusion coefficient d = D/vq. Note that 
D has also the dimension of a frequency so that d is actually dimensionless. We also 
introduce a scaled observation kernel K 1 such that -fT(:ro|a;'|) — K'{\x'\). Typically, if K 
is the indicator function of the ball of radius R, K' is the indicator of the ball of radius 
R' = R/xq. The assumption that the interaction is non local means that R 1 = 0(1). In 
other words, the observation radius is of the same order as the mean distance travelled 
by the particles between two interactions. This appears consistent with the behaviour of 
a fish, but would probably require more justifications. In the present work, we shall take 
this fact for granted. 

Let now t' = t/t , x' = x/xq the associated dimensionless time and space variables. 
Then, system (13.111) . coupled with (13.91) . (13.101) is written in this new system of units 
(after dropping the primes for the sake of clarity): 



The system now depends on only one dimensionless parameter d and two dimensionless 
functions which describe the behaviour of the fish: ^(cos^) and K(x), which are all 
supposed to be of order 1. 

Up to now, the system has been written at the microscopic level, i.e. at time and 
length scales which are characteristic of the dynamics of the individual particles. Our 




(4.1) 
(4.2) 



(4.3) 
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goal is now to investigate the dynamics of the system at large time and length scales 
compared with the scales of the individuals. For this purpose, we adopt new time and 
space units t = t /e, x = x /e with e <C 1. Then, a set of new dimensionless variables 
is introduced x = ex, t = et. In this new set of variables, the system is written (again, 
dropping the tildes for clarity): 



e{d t f e + u ■ VJ £ ) = -V, ■ (F £ f £ ) + dAuf*, 
F £ (x, lu, t) = v(uo ■ Cu £ ) (Id — to ® lu)uj £ (x, OJ, t) 

p , < J ( X . t J 

U! [X, CO, t) 



, . J £ (x,t)= / K 



x-y 



e 



(4.4) 
(4.5) 

v f{y,v,t)dydv , (4.6) 



Our goal in this paper is to investigate the formal limit e — >• of this problem. 

Our first task, performed in the following lemma, is to provide an expansion of Co 6 in 
terms of e. 



Lemma 4.1 We have the expansion: 

cu £ (x,uj,t) = tt £ (x,t) + 0(e 2 ) 

where 

Q°(x,t)= P^r, , and f(x,t) = 



\3 e (x,t)[ 



v f £ (x, v, t) dv . 



(4.7) 



(4.i 



vet 



The proof of this lemma is elementary, and is omitted. That the remainder in (14.71) is 
of order e 2 is linked with the fact that the observation kernel is isotropic. If an anisotropic 
kernel had been chosen, such as one favouring observations in the forward direction, then 
a term of order e would have been obtained. This additional term would substantially 
change the dynamics. We leave this point to future work. 

The quantity j £ (x, t) is the particle flux. We will also use the density, which is defined 
as a moment of / as well: 

p £ (x,t)= [ f £ (x,v,t)dv. (4.9) 

Thanks to lemma H~H system fl4.4p - fl4.6l) is written 

e(d t f + u ■ V x f) = -V. ■ (F £ f £ ) + dA u f + 0(e 2 ), (4.10) 
F £ (x, u, t) = v{uj ■ Q £ ) (Id - u ® u)Q £ (x, t), (4.11) 

*) = rrr% > and ^'^ *) = / v v > *) dv ■ ( 4 - 12 ) 

We note that observing the system at large scales makes the interaction local and that 
this interaction tends to align the particle velocity to the direction of the local particle 
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flux. This interaction term is balanced at leading order by the diffusion term which tends 
to spread the particles isotropically on the sphere. Obviously, an equilibrium distribution 
results from the balance of these two antogonist phenomena. 

In the remainder of the paper, we write F[f £ ] for Ffi. We introduce the operator 

Q(f) = -V„-(F[f]f) + dA ul f, (4.13) 
F[f] = is(ld-uu®uo)n[f], (4.14) 

m = rrh, and j[f]= f ufdu. (4.15) 

We note that Q[f] is a non linear operator of /, and so are F[f] and Q(f). In the 
remainder, we will always suppose that / is as smooth and integrable as necessary. We 
leave the question of finding the appropriate functional framework to forthcoming work. 

The operator Q acts on the angle variable to only and leaves the other variables x and 
t as parameters. Therefore, it is legitimate to study the properties of Q as an operator 
acting on functions of uo only. This is the task performed in the following section. 



4.2 Properties of Q 

We begin by looking for the equilibrium solutions, i.e. the functions / which cancel Q. 
Let /j, = cos 9. We denote by a(ji) an antiderivative of u(/jl), i.e. (da / dp)(n) = u(fj). We 
define 



M n (uo) =Cexp 



a(uo ■ Q) . 
d 



M n (uo)duo = 1. 



(4.16) 



The constant C is set by the normalization condition (second equality of ( 14. 16ft ) ; it 
depends only on d and on the function a but not on Q. 
We have the following: 



Lemma 4.2 (i) The operator Q can be written as 



Q(f) =dV u 



f 



(4.17) 



and we have 

H(f) :- 



gs 2 M n[f] 



duo = —d 



M 



e§ 2 



/ 



M 



duo < 0. (4.18) 



(ii) The equilibria, i.e. the functions f(oo) such that Q(f) = form a three-dimensional 
manifold £ given by 



£ = { P M n (uo) | pe 



Q e § 2 } 



(4.19) 
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and p is the total mass while Q is the director of the flux of P Mq(oj), i.e. 



[ pM n {uj)du = p (4.20) 

n = TTTTu > 3\pM a ] = [ P M n (u)udw. (4.21) 
\j[pMn\\ J we &i 

Furthermore, H(f) = if and only if f = pM^ for arbitrary p e M + and fl e § 2 . 

The function a being an increasing function of p (since v > 0), is maximal for 
• £1 = 1, i.e. for u pointing in the direction of Q. Therefore, Q plays the same role as the 
average velocity of the classical Maxwellian of gas dynamics. The role of the temperature 
is played by the normalized diffusion constant d : it measures the 'spreading' of the 
equilibrium about the average direction Q. Here the temperature is fixed by the value of 
the diffusion constant, in contrast with classical gas dynamics where the temperature is a 
thermodynamical variable whose evolution is determined by the energy balance equation. 

An elementary computation shows that the flux can be written 

j[pM n ] = (cos9) MP n, (4.22) 

where for any function ^(cos^), the symbol (g(cos9)) m denotes the average of g over the 
probability distribution Mm, i.e. 

r r n (cos 8) cxp f <t(cos e) ) dn8d8 

(g(cos8)) M = / M n (u)g(u ■ fi) duo = ' , ' . (4.23) 

Wl J/ J Km } j; e xp(^^) sin 6d6 

We note that (g (cos 0))m does not depend on Q but depends on d. In particular, 
(g(cos8)) m — > g(l) when d — >• while (g(cos8))M g, the arithmetic average of g 
over the sphere, when d —>■ oo (with g = J g(u-Q) du — | g(cos8) sin8dd). Therefore, 
(cos8)m — > 1 when cZ — ^ and (cos8)m — ► when d — > oo. For a large diffusion, the 
equilibrium is almost isotropic and the magnitude of the velocity tends to zero while for 
a small diffusion, the distribution is strongly peaked in the forward direction and the 
magnitude of the velocity tends to 1, which is the velocity of the individual particles. 

Proof of lemma l4.2t To prove (i), we introduce a reference frame such that e% = Q[f]. 
In spherical coordinates, we have 

M m (u(0,<p)) = Cexp(d-V(cos0)). (4.24) 

Therefore, 

V w (lnM n[fl ) = V w [ln{Cexp(d-V(cos#))}] 
= tTV^Hcostf)) 
= —d~ 1 v(cos8)sm9ee 

= d~ l F[f], (4.25) 
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where In denotes the logarithm and the last equality results from (13.5j) . Then, we deduce 
that 



M m V u 



f 



Mi 



dV u -[V u f-fV u QnM m , 
dA w f-V„-(F[f]f) = Q(f). (4.26) 



(14.181) follows directly from (I4.17P and Stokes theorem. 

(ii) follows directly from (i). If Q(f) = 0, then H(f) = 0. But H(f) is the integral 
of a non-negative quantity and can be zero only if this quantity is identically zero, which 
means / = pMci\f\ for a conveniently chosen p. Since ft[/] can be arbitrary, the result 
follows. The remaining statements are obvious. ■ 

Our task now is to determine the collision invariants of Q, i.e. the functions ip(uo) such 
that 

/ Q(f)*l>cLu = 0, V/. (4.27) 
Using (14.171) . this equation can be rewritten as 

/ 1 J—V u -(M m V u ,l>)du = 0, V/. (4.28) 

Jue$ 2 M n[f] 

Clearly, if ip — Constant, ip is a collisional invariant. On the other hand, there is no 
other obvious conservation relation, since momentum is not conserved by the interaction 
operator. The constants span a one-dimensional function space, while the set of equilibria 
is a three-dimensional manifold. So, we need to find some substitute to the notion of 
collisional invariant, otherwise, in the limit e — > 0, the problem will be under-determined, 
and in particular, we will lack an equation for Q (appearing in the expression of the 
equilibrium) . 

To solve the problem, we slightly change the viewpoint. We fix ft e § 2 arbitrarily, 
and we ask the problem of finding all ^'s which are collisional invariants of Q(f) for all / 
with director fl[f] = fi. Such a function ip is not a collisional invariant in the strict sense, 
because (I4.27P is valid for all / but only for a subclass of /. But this weaker concept of 
a collisional invariant is going to suffice for our purpose. So, for fixed fl, we want to find 
all ^'s such that 

/ TT y " ■ ( M ^ V -V0 du = 0, V/ such that £l[f] = ft. (4.29) 

Now, saying that fl[f] = Q is equivalent to saying that j[f] is aligned with ft[/], or again 
to 

= Q x j[f] = I /(fixw) duo. (4.30) 
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This last formula can be viewed as a linear constraint and, introducing the Lagrange 
multiplier (3 of this constraint, (3 being a vector normal to Q, we can restate the problem 
of finding the 'generalized' collisional invariants (14.291) as follows: Given Q G § 2 , find all 
such that there exist (3 G M 3 with Q ■ (3 = 0, and 

/ -^-{V UJ -(M n V^)-(3-(Qxu)M n }du = 0, V/. (4.31) 

Now, (14.311) holds for all / without constraint and immediately leads to the following 
problem for ip: 

V w ■ (MnV^V) =(3 -{fix u)M Q . (4.32) 

The problem defining ip is obviously linear, so that the set Cn of generalized collisional 
invariants associated with the vector Q is a vector space. It is convenient to introduce 
a cartesian basis (ei,e2,f2) and the associated spherical coordinates (9,(f>). Then (3 = 
(/?i,/?2,0) and (3 ■ (Q x u) = (— /?isin0 + ^ cos 0) sin 6*. Therefore, we can successively 
solve for ipi and ip2, the solutions of (14.321) with right-hand sides respectively equal to 
— sin sin 6Mq and cos 4> sin 8Mq . 

We are naturally looking for solutions in an L 2 (S 2 ) framework, since ip is aimed at 
constructing marcroscopic quantities by integration against / with respect to uj. There- 
fore, one possible framework is to look for both / and ip in L 2 (S 2 ) to give a meaning to 
these macroscopic quantities. We state the following lemma: 

Lemma 4.3 Let x £ £ 2 (§> 2 ) such that j xdu = 0. The problem 

V w • (M n V^) = X, (4.33) 

o 

has a unique weak solution in the space if 1 (§ 2 ) ; the quotient of the space if 1 (S 2 ) by the 
space spanned by the constant functions, endowed with the quotient norm. 

Proof: We apply the Lax-Milgram theorem to the following variational formulation of 
(KM): 

/ MaV u if>-V u (pdu= / x<pdu, (4.34) 

o 

for all (p G H 1 (E> 2 ). The function Mq is bounded from above and below on § 2 , so the 

o 

bilinear form at the left-hand side is continous on H 1 (E> 2 ). The fact that the average of 

o 

X over § 2 is zero ensures that the right-hand side is a continuous linear form on iJ x (S 2 ). 
The coercivity of the bilinear form is a consequence of the Poincare inequality: 3C > 

such that G H 1 ^ 2 ): 

IW > CM\o 2 := Cmm\\ij + K\\ L 2, (4.35) 
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where \ip\m is the H 1 semi-norm. We note that the Poincare inequality would not hold 
without taking the quotient. ■ 



So, to each of the right-hand sides x = ~ sin0sin$MQ or x — cos sin ^M^ which 
have zero average on the sphere, there exist solutions ipi and ip2 respectively (unique up 
to constants) of problem (14. 33D . We single out unique solutions by requesting that ipi an d 
ip2 have zero average on the sphere: J ipkdw = 0, k = 1,2. We can state the following 
corollary to lemma I4.3t 

Proposition 4.4 The set Cq of generalized collisional invariants associated with the vec- 
tor Q which belong to H 1 ^ 2 ) is a three dimensional vector space Cq = Span{l, ^2} 

More explicit forms for ipi and ip2 can be found. By expanding in Fourier series with 
respect to (ft, we easily see that 

■01 = — g(cos9) sin0, ip 2 = g(cos9) cos0, (4.36) 

where g(n) is the unique solution of the elliptic problem on [—1, 1]: 

-(1 - fi 2 )dp(e°W d (l - fi^g) + e aM/d g = -(1 - /j 2 ) 3/2 e a ^ /d . (4.37) 

We note that no boundary condition is needed to specify g uniquely since the operator 
at the left-hand side of (14.371) is degenerate at the boundaries /i = ±1. Indeed, it is an 
easy matter, using again Lax-Milgram theorem, to prove that problem (14.371) has a unique 
solution in the weighted H 1 space V defined by 

V = {g I (1 - ti 2 )- 1/2 g e L 2 (-l, 1), (1 - ^ 2 d,g e L 2 (-l, 1)}. 

Furthermore, the Maximum Principle shows that g is non-positive. 

For convenience, we introduce h{n) = (1 — ji 2 )~ l l 2 g E L 2 {— 1,1) or equivalently 
h(cos9) = g(cos9)/ sin 9. We then define 

•${u) = {VLxu) h(Q ■ u) = foe! + ip 2 e 2 ■ (4.38) 

ijj is the vector generalized collisional invariant associated with the direction Q. 

4.3 Limit e -> 

The goal of this section is to prove theorem ll.lt 

Again, we suppose that all functions are as regular as needed and that all convergences 
are as strong as needed. The rigorous proof of this convergence result is outside the scope 
of this article. 

We start with eq. (I4.10p which can be written 

e(d t f E + u ■ V x f) = Q(f) + 0(e 2 ). (4.39) 
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We suppose that f e — > / when e — > 0. Then, from the previous equation, Q(f e ) = 0(e) 
and we deduce that Q(f) = 0. By lemma fl~2l / = pMo, with p > and f2 e S 2 . Now, 
since Q operates on the variable only, this limit does not specify the dependence of / 
on (x,t), and consequently, p and Q are functions of (x,t). 

To find this dependence, we use the generalized collisional invariants. First, we con- 
sider the constant collisional invariants, which merely means that we integrate ( 14.391) with 
respect to uj. We find the continuity equation 

d t p £ + V x -f = 0, (4.40) 

where p e and j e are the density and flux as defined above. It is an easy matter to realize 
that the right-hand side is exactly zero (and not 0(e 2 )). In the limit e — > 0, p e — > p and 
j e — > j = cipQ with 

Cl = (cos^) M , (4.41) 

and we get 

d t p + V, • (cipfi) = 0. (4.42) 



Now, we multiply (14.391) by ip £ = h(u ■ Q[f £ }) (Q[f £ ] x w), integrate with respect to 
uj and take the limit e — > 0. We note that fi[/ e ] — > and that ip e is smooth enough 
(given the functional spaces used for the existence theory), and consequently, ip e — > ip = 
h(u ■ Q) (fl x u). Therefore, in the limit e — > 0, we get: 

QxX = 0, X := j (d t (pM n )+uj-V x (pM n ))h(uJ-n)ujduj. (4.43) 

Saying that Q x X = is equivalent to saying that the projection of X onto the plane 
normal to Q vanishes or in other words, that 

(Id - n ® Q)X = . (4.44) 

This is the equation that we need to make explicit in order to find the evolution equation 
for a 

Elementary differential geometry gives the derivative of with respect to Q acting 
on a tangent vector dQ to the sphere as follows: 



We deduce that 



(dn) = dr x v(u} ■ fi) (u ■ dn) M n . (4.45) 



d t (pMn) = M n (d tP + d- l vp(uj ■ d t n)), (4.46) 
(uj ■ V x )(pM n ) = M n ((uj ■ V x )p + d- x vpuj ■ ((uj ■ V x )n)). (4.47) 
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Combining these two identities, we get: 

d t ( P M n )+uj-V x (pM n ) = 

= M n [d t p + uj ■ V x p + d~ 1 vp{ u ■ d t Q + (uj ® uj) : V z ft ) ] , (4.48) 

where the symbol ':' denotes the contracted product of two tensors (if A = (Aij)ij = i 7 _ t s 
and B = (-8^)1^=1,. ..,3 are two tensors, then A : B — J2ij=i 3 Aj-By) an d V^ft is the 
gradient tensor of the vector ft: (V x Q)ij = d x .Qj . Therefore, the vector X, is given by: 

X = [dtp + u- V x p + d~ l vp(uj ■ d t Vt + (uj®oj) : V^ft)] ujhM n duj (4.49) 

The four terms in this formula, denoted by X\ to X 4 , are computed successively using 
spherical coordinates (9, <fi) associated with a cartesian basis (ei, ft) where e\ and e 2 are 
two vectors normal to ft. In the integral (14.491) . the functions ft = ft(cos#), v = u(cos9) 
and Mq = C exp ( CT(c ° s ^ ) only depend on 9. Therefore, the integrals with respect to <fi 
only concern the repeated tensor products of uj. 
We first have that f Q oj d(j) = 2tt cos 9 ft, so that 

Xl= f dtpujhMnduj = 2nd t p [ cos9h(cos9) M n (cos9) sin9d9tt, (4.50) 
Jwes 2 Jo 

and (Id - ft ® n)X 1 = 0. 

Now, an easy computation shows that 

p2tt 

/ uj g) uj d<p = 7r sin # (Id -ft (gift) + 27rcos 2 flft ®ft. (4.51) 
Jo 

We deduce that 

X 2 = {{uj®uj)V x p)hM Q duj = 
Joj& 2 

= 7r / sin 2 ft M n sin 6 1 d9 (Id - ft ® ft) V x p + 

/*7T 

+ 2vr / cos 2 ft sin 9d9 (ft ■ V^p) ft, (4.52) 

which leads to: 

(Id - ft ® ft)X 2 = n sin 2 ft M n sin d# (Id - ft <g> ft) V x p, (4.53) 
Jo 

Using (14.511) again, we find: 

X 3 = dT 1 p I ((uj®uj)d t Q)vhM n duj = 

Juj& 2 

= nd- 1 p / sin 2 9 vhM n sin9d9 (Id - ft ® ft)<9 t ft + 

/*7T 

+ 2nd~ 1 p / cos 2 v ft M c sin d9 (ft ■ d t ft) ft. (4.54) 
Jo 
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The second term at the r.h.s. of (14.541) vanishes since dt^l is normal to Q (Q being a unit 
vector). For the same reason, (Id — Q <g> Q)dtQ = dt^l and we are left with: 

/*7T 

(Id - fl ® tt)X 3 = ixd- l p I sin 2 9vhM n sm9 d6 d t VL. (4.55) 

Jo 

We now need to compute the integral with respect to (f) of the third tensor power of 
uj. After some computations, we are left with 



2?r 

2 







uj(g)uj(g)Lud(f) = iT sin 2 9 cos 9 ((Id - il <g> fi) ® Q + Q ® (Id - ® fi) + 

+ [(Id - O <g> ft) <g> ft ® (Id - ft ® n)] :2 4) 
+2tt cos 3 9 ft ® ft ® ft, (4.56) 

where the index ': 24' indicates contraction of the indices 2 and 4. In other words, the 
tensor element (J^ u <8> a; <S> ud4>)ijk equals 7r sin 2 9 cos # when (ij,k) equals any of the 
triples (1,1,3), (2,2,3), (3,1,1), (3,2,2), (1,3,1), (2, 3, 2), equals 2vr cos 3 9 when (ij, k) = 
(3, 3, 3) and is equal to otherwise. Using Einstein's summation convention, the following 
formula follows: 

( / uj ® uj ® uj d(f))V ' X VL = I / (jJ ® uj <g> uj dcf) I d Xj Qk — 
Jo \Jo J ijk 3 

= it sin 2 9 cos 9 ((Id — ft <g> ft^-ftf-S^ft/. + Sl^ (Id — ft <g> VL)j k d Xj Vt k + 

+(id-o®Q) ife %a^ fe ) 

+27rcos 3 6' fij^-fifc^fifc, (4.57) 

But since ft is a unit vector, Q k d Xj Q k = \d Xj (\Q\ 2 ) = and the first and fourth terms in 
the sum vanish. The expression simplifies into: 

/•27T 

( / uj <g> uj <g> ud(f))V x Cl = tc sin 2 # cos ((Id- ft® ft) : (V x ft))ft + 

+7r sin 2 cos# (Id - ft ® ft)((ft • V)ft), (4.58) 

The first term is parallel to ft. Besides, since ft is a unit vector, (ft • V)ft is normal to ft. 
So, we finally get 

/■27T 

(Id-ft®ft)((/ Va.fi) = vr sin 2 6» cos 6» (ft ■ V)ft, (4.59) 

Jo 

This leads to the following formula for X 4 : 

(Id-ft®ft)Xt = d _1 p (Id - ft ® ft) ( / (a; ® a; ® a;) ( V^ft) f h Mn dw J 

/"7T 

= Ttd^p sin 2 9cos9uhM n sin9d9 (ft ■ V)ft (4.60) 
Jo 



Now, we insert the expressions of X\ to X4 into (14.44ft . Using notation (14.231) . we 
finally find the evolution equation for Q: 

<T x p (sin 2 9 v h) M d t VL + cT 1 p(sin 2 6cos6v h) M (Vt ■ V)fi + 

+ (sin 2 9 h) M (Id - tt <g> ft) V x p = 0. (4.61) 

By the maximum principle, the function h is non-positive. Therefore, we can define 
similar averages as ( 14.231) . substituting Mq with sin 2 9 v h and we denote such averages 
as (<?) (sin 2 6)vhM ■ With such a notation, ( 14.611) becomes: 

p (d t VL + c 2 (tt ■ V)0) + A (Id — O ® n)V x p = 0, (4.62) 

with 

c 2 = (cos 0) (sin 2 e)i , hM , A = d ( - \ (4.63) 

\ V I (sin 2 6)vhM 

Collecting the mass and momentum eqs (I4.42p and (14.621) . we find the final macroscopic 
model of the Couzin-Vicsek algorithm: 

d t p + V, ■ (dpfi) = 0. (4.64) 
p (d t n + c 2 {yi ■ + A (Id - tt <g> fi) V x .p = 0, (4.65) 

with the coefficients ci, C2 and A given by ( 14.411) and (14.631) . This ends the proof of 
theorem 11.11 

4.4 Hyperbolicity 

The detailed study (both theoretical and numerical) of the properties of the continuum 
model (II. 4p . (11.51) . will be the subject of future work. As a preliminary step, we look at 
the hyperbolicity of the model. 

First, thanks to a temporal rescaling, t = t'/ci, we can replace c x by 1, c 2 by c := c 2 /ci 
and A by A' = A/ci. We will omit the primes for simplicity. Then, the system reads: 

dtp + V, ■ (pQ) = 0. (4.66) 

p (d t n + c(n ■ v)fi) + a (id - n ® n)v x p = o, (4.67) 

This rescaling amounts to saying that the magnitude of the velocity of the individual 
particles is equal to 1/ci in the chosen system of units. 

We choose an arbitrary fixed cartesian coordinate system (f^i, ^2,^3) and use spherical 
coordinates (9, 0) in this system (see section |3]). Then, Q = (sin 6 cos 0, sin 9 sin 0, cos#). 
A simple algebra shows that (p, 9, <fi) satisfy the system 

dtp + d x (psin9 cos0) + ^(psin^sin^) + d z (pcos9) = 0. (4.68) 
d t 9 + c(sin 9 cos d x 9 + sin 9 sin d y 9 + cos 9d z 9) + 

+A (cos 9 cos d x In p + cos 9 sin d y In p — sin 6* <9 2 In p) = 0. (4.69) 
<9 t + c(sin 9 cos 9 X + sin 9 sin d y 4> + cos ^c^) + 

+A (— sin 9 sin In p + sin 9 cos In p) = 0. (4.70) 
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Supposing that p, 9, <p are independent of x and y amounts to looking at waves which 
propagate in the z direction at a solid angle (9, 0) with the velocity director Q. In this 
geometry, the system reads: 

d t p + cos 9 d z p - p sin 9 d z 9 = 0. (4.71) 
d t 9 + ccos9d z 9- \ sin9 d z \np = 0. (4.72) 
d t (p + ccos9d z (f) = 0. (4.73) 



This is a first order system of the form 




d z p 

A(p,9,<P) | d z 9 | =0 (4.74) 

d z (p 



with 




cos 9 —p sin 9 

A(p,9, ( f ) )= ( -±f* ccos9 ] , (4.75) 


The eigenvalues 7± and 7 of the matrix A(p, 9, 0) are readily computed and are given 

by 

7o = ccos#, 7± = - Uc+ l)cos#± ((c- l) 2 cos 2 # + 4A sin 2 fl) 1/2 j . (4.76) 
2 

Two special cases are noteworthy. The case 9 = (modulo ir) corresponds to waves 
which propagate parallel to the velocity director. In this case, two eigenvalues are equal: 
7o = 7 + = c and 7_ = 1. The eigenvectors corresponding to these three eigenvalues are 
respectively the density p, and the angles 9 and 0. So far, the relative magnitude of c 
and 1 are not known. But, whatever the situation (c bigger or smaller or even equal to 
1), the matrix is diagonalizable and therefore the system is hyperbolic. 

The case 9 = tt/2 (modulo it) corresponds to waves propagating normally to the 
velocity director. In this case, 7± = ±2a/A are opposite and 70 = 0. The system for 
(p, 9) reduces to a special form of the nonlinear wave equation. The sound speed which 
propagates in the medium due to the interactions between the particles has magnitude 
equal to 2VX 

If 9 has an arbitrary value, then, a combination of these two phenomena occurs. For 
the two waves associated with 7±, there is a net drift at velocity (c + 1) cos# and two 

sound waves with velocities ((c — l) 2 cos 2 # +4A sin 2 ^) 1 ^ 2 . However, the speed of the 
wave associated with 70, is not equal to the drift of the two sound waves. A disymmetry 
appears which is not present in the usual gas dynamics equations. The resolution of the 
Riemann problem is left to future work. 
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5 Conclusion 



In this paper, we have studied the large-scale dynamics of the Couzin-Vicsek algorithm. 
For that purpose, we have rephrased the dynamics as a time-continuous one and have 
formulated it in terms of a kinetic Fokker-Planck equation. Then, a hydrodynamic scaling 
of this kinetic equation is introduced with small parameter e and the limit when e — > is 
considered. We show that the macroscopic dynamics takes place on a three dimensional 
manifold consisting of the density and director of the mean- velocity. Using a new concept 
of generalized collision invariant, we are able to derive formally the set of equations 
satisfied by the parameters and we prove that the resulting system is hyperbolic. 

Possible future directions involve the investigation of a limited range of vision in the 
backwards direction, the computation of the order e diffusive corrections, the incorpo- 
ration of more non-locality effects in the asymptotics and finally, the accounting of the 
other types of interactions, being of repulsive or attractive type. 
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